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Hybrid approaches based on relativistic hydrodynamics and transport theory have been success¬ 
fully applied for many years for the dynamical description of heavy ion collisions at ultrarelativistic 
energies. In this work a new viscous hybrid model employing the hadron transport approach UrQMD 
for the early and late non-equilibrium stages of the reaction, and 3-1-1 dimensional viscous hydrody¬ 
namics for the hot and dense quark-gluon plasma stage is introduced. This approach includes the 
equation of motion for finite baryon number, and employs an equation of state with finite net-baryon 
density to allow for calculations in a large range of beam energies. The parameter space of the model 
is explored, and constrained by comparison with the experimental data for bulk observables from 
SPS and the phase I beam energy scan at RHIC. The favored parameter values depend on energy, 
but allow to extract the effective value of the shear viscosity coefficient over entropy density ratio 
rj/s in the fluid phase for the whole energy region under investigation. The estimated value of 77 /s 
increases with decreasing collision energy, which may indicate that rj/s ol the quark-gluon plasma 
depends on baryochemical potential yLB- 


I. INTRODUCTION 

Ultra-relativistic heavy ion collisions allow to investi¬ 
gate the properties of strongly interacting matter under 
extreme conditions. At high temperatures and/or high 
net-baryon densities a new state of matter, the so-called 
quark-gluon plasma, QGP, is formed. The two main goals 
of heavy ion research are the exploration of the phase di¬ 
agram of quantum chromodynamics and the determina¬ 
tion of transport coefficients of this new state of matter. 

The studies of high energy heavy-ion collisions at the 
Large Hadron Collider (LHC) at CERN and the Rela¬ 
tivistic Heavy Ion Collider (RHIC) at Brookhaven Na¬ 
tional Laboratory have revealed that the quark-gluon 
plasma behaves like an almost perfect fluid. In recent 
years, so-called hybrid approaches based on (vis¬ 
cous) relativistic hydrodynamics for the hot and dense 
stage coupled to hadron transport approaches for the de¬ 
coupling stage of the reaction have been applied with 
great success to extract average values of the shear vis¬ 
cosity over entropy ratio ry/s. The results are very close 
to the conjectured universal limit of 77 /s = based 
on the anti-de Sitter/conformal field theory (AdS/CFT) 
correspondence [S]. For example, the values extracted 
in Ref. [7] are rj/s = 0.12 for collisions at RHIC, and 
77 /s = 0.2 at the LHC. 

One expects the formation of partonic matter in 
heavy ion collisions at ultra-relativistic energies (see, e.g.. 
Ref. 0). However, it is unknown at what collision energy 
the transition from hadronic to partonic matter sets in. 
In addition, as the collisions at lower energies probe the 
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phase diagram at larger net-baryon densities, it may be 
possible to experimentally see signs of the theoretically 
predicted critical point [5] and the first-order phase tran¬ 
sition beyond it. To investigate these questions the so- 
called beam energy scan (BES) programs at SPS (NA49, 
NA61 experiments) and at RHIC (STAR, PHENIX ex¬ 
periments) were started. One of the surprises of the stage 
I of the BES program at RHIC has been that the px- 
differential elliptic flow, V 2 {pt), of charged hadrons does 
not change significantly when the collision energy is re¬ 
duced from ^snn = 200 to ~ 20 CeV (TU]. The large 
values of elliptic flow measured at ^snn = 200 CeV col¬ 
lisions were taken as a sign of very low shear viscosity 
of the matter formed in these collisions. Thus, the large 
V 2 {j>t) measured in collisions at lower energy leads to 
the question how 77 /s changes as function of net-baryon 
density and baryochemical potential pb m- 

Unfortunately, many of the hydrodynamical and hy¬ 
brid models used to model collisions at full RHIC and 
LHC energies are not directly applicable to collisions at 
RHIC BES and CERN SPS energies, nor to collisions 
at even lower energies in the future at Facility for An¬ 
tiproton and Ion Research (FAIR), Nuclotron-based Ion 
Collider Facility (NICA) and the stage H of the BES 
program at RHIC. The simplifying approximations of 
boost invariance and zero net-baryon density are not 
valid, and different kinds of non-equilibrium effects play a 
larger role. To overcome these limitations, a novel hybrid 
approach has been developed. This approach is based 
on the Ultra-relativistic Quantum Molecular Dynamics 
(UrQMD) transport |I2j for the non-equilibrium early 
and late stages, and on a (3-|-l)-dimensional viscous hy¬ 
drodynamical model m for the hot and dense stage of 
the reaction. 

In this paper, this approach is applied to extract the 
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shear viscosity coefRcient over entropy density ratio of 
strongly interacting matter from the heavy-ion collision 
data at RHIC beam energy scan energies in the broad 
range ^snn = 7.7-200 GeV. The details of the model 
are explained in Section and the generic effects of 
finite shear viscosity on the hydrodynamical expansion 
are described in Section m The sensitivity of particle 
yields and spectra to the parameters for the initial and 
final state transitions is explored in Sectionjl^ Section[V| 
contains the main results of this work including the esti¬ 
mated values of the effective shear viscosity over entropy 
density ratio as a function of beam energy. Finally, the 
main conclusions are summarized and an outlook on fu¬ 
ture work is given in Section [VI| 

II. MODEL DESCRIPTION 

Our hybrid approach combines the UrQMD transport 
model [12] for the early and late stages of the evolu¬ 
tion with a dissipative hydrodynamical model, called 
vHLLE HU, for the hot and dense stage. The distinguish¬ 
ing features of the present model are that the fluid dy¬ 
namical expansion is solved numerically in all three spa¬ 
tial dimensions without assuming boost invariance nor 
cylindrical symmetry, the equations of motion for finite 
net-baryon and charge densities are explicitly included 
and, in contrast to the standard UrQMD hybrid approach 
(UrQMD-3.4 at urqmd.org) jS] ITTj . dissipation in the 
form of shear viscosity is included in the hydrodynamical 
stage. Different to our previous studies [ndH] , event- 
by-event fluctuations are now included. The hadronic 
cascade operates with the full phase-space distribution 
of the final particles which allows for a proper compari¬ 
son to experimental data. 


A. Pre-thermal Phase 

The UrQMD string/hadronic cascade is used to de¬ 
scribe the primary collisions of the nucleons, and to cre¬ 
ate the initial state of the hydrodynamical evolution. The 
two nuclei are initialized according to Woods-Saxon dis¬ 
tributions and the initial binary interactions proceed via 
string or resonance excitations, the former process be¬ 
ing dominant in ultrarelativistic collisions (including the 
BES collision energies). All the strings are fragmented 
into hadrons before the transition to fluid phase (fluidiza¬ 
tion) takes place, although not all hadrons are yet fully 
formed at that time, i.e., they do not yet have their free- 
particle scattering cross sections, and thus do not yet 
interact at all. The hadrons before conversion to fluid 
should not be considered physical hadrons, but rather 
marker particles to describe the flow of energy, momen¬ 
tum and conserved charges during the pre-equilibrium 
evolution of the system. The use of UrQMD to ini¬ 
tialise the system allows us to describe some of the pre¬ 
equilibrium dynamics and dynamically generates event- 



FIG. 1: The earliest possible starting time of the hydrody¬ 
namic evolution as a function of according to Eq. [^ 


by-event fluctuating initial states for hydrodynamical 
evolution. 

The interactions in the pre-equilibrium UrQMD evolu¬ 
tion are allowed until a hypersurface of constant Bjorken 
proper time tq = \/t'^ — is reached, since the hydro¬ 
dynamical code is constructed using the Milne coordi¬ 
nates {T,x,y,r]), where r = \/t‘^ — [13]. The UrQMD 

evolution, however, proceeds in Cartesian coordinates 
(t, X, y, z), and thus evolving the particle distributions to 
constant r means evolving the system until large enough 
time ti in such a way that the collisional processes and 
decays are only allowed in the domain — z"^ < tq. 
The resulting particles on t = ti surface are then prop¬ 
agated backwards in time to the r = tq surface along 
straight trajectories to obtain an initial state for the hy¬ 
drodynamic evolution. 

The lower limit for the starting time of the hydrody¬ 
namic evolution depends on the collision energy accord¬ 
ing to 

To = 2R/ (-^fiNN/SmAr)^ — 1, (1) 

which corresponds to the average time, when two nuclei 
have passed through each other, i.e., all primary nucleon- 
nucleon collisions have happened. This is the earliest pos¬ 
sible moment in time, where approximate local equilib¬ 
rium can be assumed. The tq values calculated according 
to this formula are plotted in Fig. [^ 

To perform event-by-event hydrodynamics using fluc¬ 
tuating initial conditions every individual UrQMD event 
is converted to an initial state profile. As mentioned, the 
hadron transport does not lead to an initial state in full 
local equilibrium, and the thermalization of the system 
at r = To has to be artihcially enforced. The energy and 
momentum of each UrQMD particle at tq is distributed 
to the hydrodynamic cells ijk assuming Gaussian density 
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FIG. 2: An example of a fluctuating (single event) initial en¬ 
ergy density profile in the transverse plane at ?7 = 0. The pro¬ 
file is obtained with R± = = 1 fm Gaussian smearing and 

corresponds to a 20-30% Au-Au collision at y^s^N = 39 GeV. 


profiles 
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where Ax^, Ay^, Arji^ are the differences between parti¬ 
cle’s position and the coordinates of the hydrodynamic 
cell fc}, and 7 ^ = cosh(?/p — rj) is the longitudinal 
Lorentz factor of the particle as seen in a frame moving 
with the rapidity rj. The normalization constant C is cal¬ 
culated from the condition that the discrete sum of the 
values of the Gaussian in all neighboring cells equals one. 
The resulting AP“ and AN^ are transformed into Milne 
coordinates and added to the energy, momentum and 
baryon number in each cell. This procedure ensures that 
in the initial transition from transport to hydrodynamics 
energy, momentum and baryon number are conserved. 

For the present study energy and momentum of the 
initial particles are converted at tq into a perfectly equi¬ 
librated fluid, i.e., the initial values for the viscous 
terms in the energy-momentum tensor are set to zero: 

= n(ro) = 0. In other words the components 
of the energy-momentum tensor stay the same, but the 
components change, when we switch from UrQMD 
to the fluid. Thus, we do not consider how much the 
energy-momentum tensor of UrQMD deviates from the 
ideal fluid energy-momentum tensor, but leave this topic 
for further studies. 

One typical example of the initial energy density dis¬ 
tributions in the transverse plane at midrapidity for one 
event is presented in Fig. The parameters R± and P,j 
regulate the granularity or the initial state. At the same 


time they influence the initial entropy of the hydrody¬ 
namic evolution, while the total initial energy and mo¬ 
mentum are always fixed to be equal to the energy and 
momentum of the pre-equilibrium UrQMD event. The 
dependence of the final results on these two parameters 
is discussed later in Section 


B. Hydrodynamic Evolution 

The (3-|-l)-dimensional viscous hydrodynamical code 
vHLLE is described in full detail in Ref. m- We repeat 
here only its main features. The code solves the usual 
local energy-momentum conservation equations 

= 0, (4) 

= 0, (5) 

where iVg and Nq are net baryon and electric charge 
currents respectively, and the semicolon denotes the co¬ 
variant derivative. The calculation^ is done in Milne 
coordinates (r,a;,y,ry), where r = \/t^ — and ry = 
l/21n[(t -I- z)/{t — z)]. 

In the Israel-Stewart framework of causal dissipative 
hydrodynamics m. the dissipative currents are inde¬ 
pendent variables. For the purpose of the present work 
we set the bulk viscosity to zero. C/s = 0. We work 
in the Landau frame, where the energy diffusion flow is 
zero, and neglect the baryon and charge diffusion cur¬ 
rents, which is equivalent to zero heat conductivity. For 
the shear stress evolution we choose the relaxation time 
Ttt = 5r]/{Ts), the coefficient = 4/3r^, and approx¬ 
imate all the other coefficients EOlllI] by zero. For the 
shear-stress tensor we obtain the evolution equation 

= -( 6 ) 

Ttt O 

where the brackets denote the traceless and orthogonal 
to part of the tensor and n/// is the Navier-Stokes 
value of the shear-stress tensor. 

Another necessary ingredient for the hydrodynamic 
stage is the equation of state (EoS) of the medium. 
We use the chiral model EoS [22], which features cor¬ 
rect asymptotic degrees of freedom, i.e., quarks and glu¬ 
ons in the high temperature and hadrons in the low- 
temperature limits, crossover-type transition between 
confined and deconfined matter for all values oi y-B and 
qualitatively agrees with lattice QCD data at yrs = 0. 

The tests to confirm the accuracy of the code have 
been reported in Ref. m- In particular the solutions 
have been compared to the ideal Gubser solution [23] 
and to a numerical solution of dissipative hydrodynamics 
calculated using the VISH2+1 hydro code [2i] . 


^ Typical grid spacing used in the calculations: Ax = Ay = 0.2 fm, 
A 77 = 0.05 — 0.15 and timestep At = 0.05 — 0.1 fm/c depending 
on the collision energy. A finer grid with Ax = Ay = 0.125 fm 
was taken to simulate peripheral collisions. 
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C. Particlization and Hadronic Rescattering 

It is well known that hydrodynamics loses its valid¬ 
ity when the system becomes dilute. To deal with this 
problem we apply the conventional Cooper-Frye prescrip¬ 
tion [25] to particlize the system (convert the fluid to 
individual particles) at a hypersurface of constant lo¬ 
cal rest frame energy density, and use the UrQMD cas¬ 
cade to describe the further evolution of these particles. 
This switching hypersurface is evaluated during the hy¬ 
drodynamic evolution using the Cornelius routine |26j . 
and as a default value for the switching density we use 
Esw = 0.5 GeV/fm^, which in the chiral model EoS cor¬ 
responds to T « 175 MeV at = 0. At this energy 
density the crossover transition is firmly on the hadronic 
side, but the density is still a little higher than the chem¬ 
ical freeze-out energy density suggested by the thermal 
models m- Thus the hadronic transport can take care 
of both chemical and kinetic decoupling of hadrons. We 
discuss the sensitivity of the results to the value of the 
switching density in section |IV[ 

As given by the Cooper-Frye prescription, the hadron 
distribution on each point of the hypersurface is 


in the 3-1-1 dimensional numerical solution of the fluid- 
dynamical equations using Milne coordinates is slightly 
violated as discussed in Refs. naiiii. 

To take into account the dissipative corrections to the 
distribution function /, we use the well-known Grad’s 
14-moment ansatz for a single component system, and 
assume that the correction is the same for all hadron 
species. We evaluate the particle distribution in the 
rest frame of the fluid at each surface element using the 
Gooper-Frye formula 


d^AN, 


dp*d{cos 9)d4> 







isotropic 


1 + (1 T /eq) 


2T2(e-fp) ■ 


( 8 ) 


Wvi 


The distribution function in Eq. § is expressed in 
terms of temperature and chemical potential(s), which 
implies a grand canonical ensemble. This allows to do 
the particle sampling independently on each surface el¬ 
ement. To create an ensemble for particles, we perform 
the following steps at each element Aai'. 


• Eirst, the average number of hadrons of every sort 
is calculated: 


The phase space distribution function / is usually as¬ 
sumed to be the one corresponding to a noninteracting 
hadron resonance gas in or close to the local thermal equi¬ 
librium. This is inconsistent with mean fields included 
in the chiral model EoS used during the evolution. To 
solve this inconsistency we evaluate the switching sur¬ 
face using the chiral model EoS, but use a free hadron 
resonance gas EoS to recalculate the energy density, pres¬ 
sure, flow velocity temperature, and chemical poten¬ 
tials from the ideal part of the energy-momentum tensor 
and charge currents, and use these values to evaluate 
the particle distributions on the switching surface. For 
example the above mentioned temperature of T « 175 
MeV in chiral model EoS at zero baryon density and 
Csw = 0.5 GeV/fm^, drops to T « 165 MeV in the free 
hadron resonance gas EoS. This procedure ensures that 
the total energy of the produced particles is reasonably 
close to the overall energy flow through the particliza¬ 
tion hypersurface (up to negative contributions to the 
Gooper-Frye formula), although a small error arises since 
we use a different energy density to evaluate the posi¬ 
tion of the surface, and the properties of the fluid on 
it^. We have checked that in a case of event-averaged 
initialization, this error is on the level of few percents. 
In addition, the conservation of energy and momentum 


^ The exact procedure suggested in Ref. m requires a numeri¬ 
cal solution of a cubic equation for each surface element, and is 
therefore too slow for event-by-event studies. 


ANi = Acr^w'^rii^th = AfJoni^th 

• For a given (Wot) = the number of particles 

to be created is generated according to a Poisson 
distribution with a mean value (Wot)- 

• For each created particle, the type is randomly cho¬ 
sen based on the probabilities W/Wot- 

• A momentum is assigned to the particle in two 
steps: 

1. The modulus of the momentum is sampled in 
the local rest frame of the fluid, according to 
the isotropic part of Eq. ^ , and the direction 
of momentum is picked randomly in 47r solid 
angle. 

2. The correction for IWesidual IW-esidual ’ lAyisc 
in Eq. Q is applied via rejection sampling: 
A random number x in the range [0, Wmax] 
is generated. If a; < W, the generated mo¬ 
mentum is accepted, if not, the momentum 
generating procedure is repeated. 

• The particle momentum is Lorentz boosted to the 
center of mass frame of the system. 

• The particle position is taken to be equal to the co¬ 
ordinate of the centroid of the corresponding sur¬ 
face element, except for the spacetime rapidity of 
the particle, which is uniformly distributed within 
the longitudinal size of the volume element. 
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FIG. 3: (Color online) Transverse momentum spectra of neg¬ 
ative pions, positive and negative kaons in Euh = 40 A GeV 
(■\/snn = 8.8 GeV) central Pb-Pb collisions. The experimen¬ 
tal data from the NA49 collaboration |29| are compared to 
the hybrid model calculations with rj/s = 0 (dashed lines) 
and rj/s = 0.2 (solid lines) in the hydrodynamic phase. The 
results from UrQMD model with no intermediate hydro phase 
(dubbed as “pure UrQMD”) are shown with dotted lines. 


For the current study, no correction over the grand 
canonical procedure is made to effectively account for the 
exact conservation of the total baryon/electric charge, 
strangeness and total energy/momentum for every sam¬ 
pled event^. As a result, these quantities fluctuate from 
event to event. 

The generated hadrons are then fed into the UrQMD 
cascade. Since the cascade accepts only a list of parti¬ 
cles at an equal Cartesian time as an input, the created 
particles are propagated backwards in time to the time 
when the first particle was created. The particles are not 
allowed to interact in the cascade until their trajectories 
cross the particlization hypersurface. 


III. SENSITIVITY TO SHEAR VISCOSITY 

The overall effects of shear viscosity on hydrodynami- 
cal expansion have been extensively discussed in the lit¬ 
erature [231 . Here we show that the results from 

high energy collisions, e.g., entropy increase, enhance¬ 
ment of transverse and inhibition of longitudinal expan¬ 
sion, and suppression of anisotropies are also manifested 
in calculations at lower collision energies. 

We have carried out event-by-event simulations for dif¬ 
ferent collision energies, centralities, and two fixed values 
of shear viscosity: rj/s = 0 (ideal hydro evolution) and 
77 /s = 0.2. For these simulations we use the values of the 


® For a suggested procedure to impose the conservation laws, see 
Ref. [26]. 




FIG. 4: (Golor online) Pion and kaon dN/dy in Eiab = 
40 A GeV (Usnn = 8.8 GeV) central Pb-Pb collisions (top) 
and charged hadron dN/drj distributions at ^snn = 19.6, 39, 
62.4 and 200 GeV central Au-Au collisions (bottom). The ex¬ 
perimental data from the NA49 [29] and the PHOBOS collab¬ 
orations m are compared to the hybrid model calculations 
with 77 /s = 0 (dashed lines) and rj/s = 0.2 (solid lines) in the 
hydrodynamic phase. 

Gaussian radii for the particles’ energy/momentum depo¬ 
sition R± = = I fm (see Eqs. ® and ^). The initial 

time is chosen according to Eq. (1H, however for the col¬ 
lisions at energies equal or higher than y'sNN = 27 GeV 
we set To = 1 fm/c. 

To reduce the need for CPU time, we use so called 
oversampling technique, as in Ref. [33]. For each colli¬ 
sion energy, centrality, and parameter set we have created 
around 500 hydrodynamic events with randomly gener¬ 
ated initial conditions. For each hydrodynamic event, 
or transition hypersurface, we generate Aoversampie = 
50— 100 final state events, which results in 25000 — 50000 
events used to calculate observables. We have checked 
that the oversampling procedure does not significantly 
affect the final observables by creating 1000 or 10000 
hydrodynamic events, with A"oversampie = 20 and 2 re¬ 
spectively, for several parameter sets. In both cases the 
calculated observables agreed within statistical errors. 

The available experimental data set for the basic bulk 
hadron observables at the BES energies is inhomoge- 
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FIG. 5: (Color online) pt integrated (0.2 < pr < 2.0 GeV and 
\ri\ < 1 ) elliptic ( 112 ) and triangular (vs) flows of all charged 
hadrons in 20-30% central Au-Au collisions as a function of 
collision energy, calculated with the event-plane method. The 
elliptic and triangular flow data is from the STAR collabo¬ 
ration [M SO]. The solid line depicts the calculation with 
p/s = 0 .2, the dashed line with 77 /s = 0 whereas the dotted 
line corresponds to the “pure” UrQMD calculation with no 
intermediate hydrodynamic stage. 



neons. (Pseudo)rapidity spectra of all charged hadrons 
for Au-Au collisions are available from the PHOBOS 
analysis [31] for ^snn = 19.6, 62.4 and 200 GeV en¬ 
ergies only. The pp spectra are published for y^SNN = 
62.4 GeV by the PHOBOS collaboration [35] and for 
= 200 GeV by the PHENIX collaboration [53] . 
To cover the lower collision energies we use dN/dy and 
PT-spectra from the NA49 [35] collaboration for Pb-Pb 
collisions at Euh = 40 and 158 A GeV, which correspond 
to y'SNN = 8.8 and 17.6 GeV, and set up the simula¬ 
tions accordingly for Pb-Pb system. For the elliptic flow 
we compare to the STAR results at ^snn = 7.7, 11.5, 
19.6, 27, 39 GeV [10] and 200 GeV [37] collision energies. 
In the model we define the centrality classes as impact 
parameter intervals based on the optical Glauber model 
estimates [38l[39] . 

The transverse momentum distributions of identified 
particles at ^snn = 8.8 GeV (Aiab = 40 A GeV) colli¬ 
sion, and (pseudo) rapidity distributions of identified or 
charged particles at collision energies y'SNN = 8.8-200 
GeV are shown in Figs. [^ and [^ respectively. As can 
be seen, the inclusion of shear viscosity in the hydro- 
dynamic phase hardens the pp spectra, and increases 
dN/dy (and similarly dN/dp) at midrapidity, squeezing 
the overall rapidity distribution. This effect can be at¬ 
tributed to the effect of shear viscosity on the strong 
longitudinal expansion of the system in the initial state 
for the hydrodynamic phase. Shear viscosity attempts 
to isotropize the expansion by decelerating it in the lon¬ 
gitudinal direction and accelerating it in the transverse 
direction. The energy of the hydrodynamic system is al¬ 
ways conserved, whereas additional entropy is produced 
during the viscous hydrodynamic evolution, which ex¬ 


plains the increased total particle multiplicity. Gompar- 
ing to the experimental data one observes that p/s = 0.2 
gives a good estimate of the rapidity and transverse mo¬ 
mentum distributions at the lowest collision energy point 
^snn = 8.8 GeV (Aiab = 40 A GeV), while overestimat¬ 
ing dN/dp at midrapidity for the rest of collision energies 
except for the highest energy, ^snn = 200 GeV, where 
we underestimate the PHOBOS results. 

In Fig. [^ the pT-averaged elliptic and triangular flow 
coefficients V 2 and U 3 are shown as a function of collision 
energy. The flow coefficients are calculated using the 
event-plane method as described in Ref. [33], including 
the event plane resolution correction. As expected, the 
elliptic and triangular flow coefficients are suppressed by 
the shear viscosity. However, comparing the results for 
p/s = 0.2 to the STAR experimental results at 20 — 30% 
centrality we find that the suppression is too weak for 
y^snn 30 GeV and too strong otherwise. The latter is 
consistent with the fact that the optimal value of p/s re¬ 
quired to fit the elliptic flow data at ^snn = 200 A GeV 
is p/s = 0.08 assuming the initial energy density profile 
from Monte Carlo Glauber approach [51]. Another par¬ 
ticular feature of both V 2 curves is that, in the region 
y^SNN ~ 20-62 GeV the elliptic flow decreases as a func¬ 
tion of .^snn- If we do not limit the initial time tq from 
below at energies ^snn > 25 GeV, but take it directly 
from Eq. ([^, we do not see this decrease, but V 2 increases 
monotonously with increasing collision energy. Thus we 
expect that the reason for the nonmonotonous behavior 
is in our choice for the initial time of the hydrodynamic 
evolution. 

The results from the standard UrQMD cascade (with¬ 
out intermediate hydrodynamic phase) are also shown for 
comparison on Figs. ]^ and [^ with dotted lines. One may 
conclude that, whereas standard UrQMD does a good 
job for pT-spectra and rapidity distributions at the low¬ 
est energy, it clearly underestimates V 2 when the collision 
energy increases (which repeats the conclusion about the 
V 2 excitation function from Ref. [53], and later results 
from V 3 analysis in Ref. M)- This is an indication of 
too large viscosity of the high-density medium and served 
historically as a motivation to introduce the intermediate 
hydrodynamic stage. 


IV. INVESTIGATION OF PARAMETER SPACE 


After investigating the generic influence of a finite 
shear viscosity during the hydrodynamic evolution on ba¬ 
sic bulk observables, it is clear that we cannot fit all the 
available experimental data using the same set of param¬ 
eters'^. Thus we have to adjust the model parameters 


^ The internal parameters of UrQMD, e.g., particle properties 
and cross sections, are fixed by experimental data as explained 
in Ref. US). The effects of changes in resonance properties were 
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relative change of the parameter 



relative change of the parameter 


FIG. 6: Parameter dependence of the total yield at midra¬ 
pidity (top) and the effective temperature of pion, kaon and 
proton pt spectra (bottom) in 0-5% central An-An collisions 
at Vsnn = 19.6 GeV. 


according to the collision energy before drawing any con¬ 
clusions about the physical properties of the system. 

In this section we study systematically the sensitivities 
of the particle yield at midrapidity, which is a measure 
for the final entropy, the effective slope parameter that 
measures the strength of the transverse expansion, and 
the anisotropic flow to the main parameters of the model. 
Due to the limited space, and to emphasize the main fea¬ 
tures of the dependencies, we restrict ourselves to one 
collision energy, y'sNN = 19-6 GeV, in the middle of the 
investigated range. Since the influence of shear viscosity 
was discussed above, we now concentrate on the remain¬ 
ing parameters of the model: the two Gaussian radii R± 
and Rjj for the initial distribution of energy, momentum 


studied in Ref. m- It was found that if the changes stay within 
experimentally acceptable limits, the effects on final particle dis¬ 
tributions are small. 



relative change of the parameter 


FIG. 7: (Color online) Parameter dependence of py integrated 
elliptic flow V 2 of charged hadrons in 20-30% central Au-Au 
collisions at ^snn = 19.6 GeV. The experimental value of the 
elliptic flow is shown with a solid red line for comparison. 


and charges, the starting time for the hydro phase tq, and 
the energy density Csw when the switch to the hadronic 
cascade happens. The default case is R± = IL = 1.0 fm. 
To = 1.22 fm/c (calculated according to Eq. [l|, 77 /s = 0 
(for simplicity) and = 0.5 GeV/fm^. The depen¬ 
dencies are presented in Figs. and where each curve 
corresponds to the variation of only one of the parame¬ 
ters, while keeping the default values for the others. All 
values are normalized to their default values to allow a 
direct comparison to each other. The effective tempera¬ 
tures of the hadron spectra in the lower panel of Fig. 
are defined as the parameter of the exponential fit: 

dN _ / rriT \ 

mTdmTdy~ ly T,s) ' 

where the rriT — m range is [0.2-1] GeV for pions and 
protons and [0.05-1] GeV for kaons®. In general we do 
observe only a very weak dependence on the parameters, 
that is less than 10 % for a 10 % change in parameters. 
The observed dependencies can be summarized as: 

• Increased i?j_ smoothens the initial energy den¬ 
sity profile in the transverse plane, which leads to 
smaller gradients and less explosive transverse ex¬ 
pansion. The latter leads to a decrease of the ef¬ 
fective temperature (inverse slope) Tgg of the pr- 
spectra, see Fig. lower panel. Larger R± also 
results in decreased ellipticity and triangularity of 
an initial energy density profile, which is hydrody- 
namically translated into smaller hnal elliptic (v 2 , 
see Fig.[^ and triangular ( 7 ) 3 ) flow components. 


® Smaller my — m range for pions and protons is taken since the 
lowest my — m part of the spectrum has a different slope than 
the intermediate my — m range. 
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TABLE I: Schematical representation of the response (in¬ 
crease or decrease) of the observables to the increase of a 
particular parameter of the model. 

• In a similar manner, the increase of i?,, leads to 
shallower longitudinal gradients and weaker longi¬ 
tudinal expansion. Thus more energy remain at 
midrapidity to form stronger transverse expansion, 
which increases Tefj and V 2 - On the other hand, 
larger also results in larger initial entropy of 
the system, which considerably increases the final 
particle multiplicity, see Fig. upper panel. 

• Increased tq has two effects: 

1. It leads to a shorter lifetime of the hydrody¬ 
namic phase, as a result of longer pre-thermal 
phase. 

2. At the same time tq enters the Gaussian en¬ 
ergy/momentum smearing profile. Thus its 
increase acts opposite to the increase of i?,,. 

From the response of the observables to the increase 
of To we find that the second effect is stronger. 

• Increased Csw shortens the effective lifetime of the 
hydrodynamic phase. The shorter time to develop 
radial and elliptic flows is not fully compensated by 
the longer cascade phase, which results in the de¬ 
crease of both final Tea and final V 2 - Since the total 
entropy is conserved in the ideal hydrodynamic ex¬ 
pansion, but increases in the cascade stage, the final 
particle multiplicity increases with the increase of 

^SW 

The observed dependencies are schematically depicted 
in Table |Tj where the signs of the responses of the observ¬ 
ables to the increase of a particular model parameter are 
shown. As for the magnitudes of the response, one can 
see from the plots that by varying the parameters of the 
initialization procedure one has a nearly linear influence 
on the final dN/dy, T^s and V 2 - From Fig. one can 
see that by choosing a larger value of R± it is possible to 
approach the experimental value of V 2 with zero shear vis¬ 
cosity in the hydrodynamic phase. However, such value 
is inconsistent with the pr spectra and charged particle 
multiplicity. 

Investigating all the dependencies in detail allows us to 
choose parameter values which lead to a reasonable repro¬ 
duction of the data. These values are shown in Table El 
For reasons of simplicity we keep Cgw = 0.5 GeV/fm^ 
for all collision energies, since the other parameters pro¬ 
vide enough freedom for adjustment. Note that since 


[GeV] 

To [fm/c] 

R± [fm] 

Rr, [fm] 

y/s 

7.7 

3.2 

1.4 

0.5 

0.2 

8.8 (SPS) 

2.83 

1.4 

0.5 

0.2 

11.5 

2.1 

1.4 

0.5 

0.2 

17.3 (SPS) 

1.42 

1.4 

0.5 

0.15 

19.6 

1.22 

1.4 

0.5 

0.15 

27 

1.0 

1.2 

0.5 

0.12 

39 

0.9* 

1.0 

0.7 

0.08 

62.4 

0.7* 

1.0 

0.7 

0.08 

200 

0.4* 

1.0 

1.0 

0.08 


TABLE IL Collision energy dependence of the model param¬ 
eters chosen to reproduce the experimental data in the BES 
region and higher RHIC energies. An asterisk denotes the 
values of starting time to which are adjusted instead of being 
taken directly from Eq. 

the model requires a lot of GPU time to obtain results 
for each particular collision energy and centrality, it is 
at the moment impractical to provide y^-optimized val¬ 
ues of the model parameters and their errors. Thus the 
parameters are adjusted manually based on a visual cor¬ 
respondence to the data. A full fledged fit to the 
data is planned for the future using a model emulator, as 
suggested in Refs. [48ll45] . 


V. RESULTS FOR BULK OBSERVABLES 

Finally, let us have a look at the results for bulk ob¬ 
servables with the energy dependent parameters for the 
hydrodynamic description (see Table 0. 

The (pseudo)rapidity spectra are presented in Fig. 
One can see that whereas the parameters were adjusted 
to reproduce the total multiplicities, the resulting shapes 
of the pseudorapidity distributions are also in a reason¬ 
able agreement with the data. From the model results 
one can observe the change in shape from the single peak 
structure at y'sNN < 20 GeV to a doubly-peaked distri¬ 
bution (or from a Dromedary to a Bactrian camel shape) 
which starts to form at .^snn = 39 GeV. At higher colli¬ 
sion energies we observe a shallow dip around zero pseu¬ 
dorapidity. 

The pt spectra of pions, kaons and protons in collisions 
at y'SNN = 62.4, 17.6, and 8.8 GeV energies are shown in 
Fig. [9} In general the spectra and especially the pt slopes 
are reproduced, which indicates that both the collective 
radial flow (generated in the hydrodynamic and cascade 
stages), and thermal motion are combined in the right 
proportion. 

The elliptic and triangular flow coefficients for 20-30% 
central Au-Au collisions as a function of collision energy 
are presented in Fig.[0j As expected, the calculated val¬ 
ues of the elliptic flow follow the data closely, since this 
quantity was used to fix the parameters. In contrast to 
that, triangular flow V 3 is calculated from the same sim- 
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FIG. 8: (Color online) Pseudorapidity distributions of 
charged hadrons (top) in Au-Au collisions at .^snn = 19.6, 39, 
62.4 and 200 GeV energies, and rapidity distributions of iden¬ 
tified hadrons in Pb-Pb collisions at i?iab = 158 and 40 A GeV 
(V®nn = 17.6 and 8.8 GeV) energies (middle and bottom 
panels, respectively). The calculations were done using the 
collision energy dependent parameters listed in Table [III The 
data are from the PHOBOS [34] and the NA49 [29] collabo¬ 
rations. 


ulated events, and thus can be considered as a prediction 
of the model. We expect that the non-monotonous be¬ 
havior of Vs is an artifact of our fitting procedure, and 
more careful adjustment of the model parameters would 


FIG. 9: (Color online) pr spectra of identified hadrons in Au- 
Au collisions at y^SNN = 62.4 GeV energy (top) and in Pb-Pb 
collisions at Slab = 158 and 40 A GeV (y^SNN = 17.6 and 
8.8 GeV) energies (middle and bottom panels, respectively). 
The model calculations were carried out using the collision 
energy dependent parameters listed in Table jllj and the data 
are from the PHOBOS and NA49 collaborations [ 29113511^ . 


further smoothen the behavior of U 3 (-\/s). 

The 20-30% centrality class was chosen because the el¬ 
liptic flow signal is strongest around this centrality class. 
Also, at this centrality nonflow contributions from mini¬ 
jets, which are not included in the model, are small. The 
centrality dependence of elliptic flow at y'SNN = 39 GeV 
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FIG. 10: (Color online) Pt integrated elliptic and triangular 
flow coefficients V 2 and V 3 as a function of collision energy. 
Both the experimental and calculated coefficients were evalu¬ 
ated using the event plane method. The calculation was done 
using the collision energy dependent parameters listed in Ta¬ 
ble]^ and the data is from the STAR collaboration [13110]. 


FIG. 12: (Color online) Effective values of shear viscosity over 
entropy density rj/s used to describe the experimental data 
at different collision energies as shown in Table [n| The green 
band represent an estimate of uncertainty in rj/s resulting 
from the allowed variation of model parameters around their 
optimal values. 



centrality percentage 


FIG. 11: (Color online) pr integrated elliptic flow coefficient 
V 2 in ^snn = 39 GeV Au-Au collisions as function of central¬ 
ity. Both the experimental and calculated V 2 was evaluated 
using the event plane method. The calculation was done using 
the collision energy dependent parameters fisted in Table |TT| 
and the data is from the STAR collaboration |10 |. 


is shown in Fig. 11 The parameters are the same at 
all centralities. In peripheral collisions the model sig¬ 
nificantly undershoots the data. This is due to the 
smoothening procedure used to convert individual parti¬ 
cles to the fluid-dynamical initial state. With the present 
smearing parameters the eccentricity of the system is too 
small in peripheral collisions, where the size of the entire 
system is comparable to the smearing radius. 

The most important conclusion from the adjustment 
procedure is that reproduction of the data requires an 
effective rj/s which decreases as a function of increasing 
collision energy, see Table [H] and Fig. [I^ On Fig. [T^ one 
can also see an estimated error band around the optimal 
values of rj/s. As mentioned, a proper determination of 


the error bars would require a fit. Currently the error 
band is estimated from the variations of two parameters 
of the model ( 77 /s and Rp) which result in the same value 
of Pt integrated elliptic flow and a 5% variation in the 
slope of proton pp spectrum, which is the most sensitive 
to a change in radial flow. 

In the present calculations 77 /s is taken to be con¬ 
stant during the evolution of the system, and its value 
changes only with the collision energy. However we ex¬ 
pect that physical 77 /s depends on both the temperature 
and baryon chemical potential, and that p/s has a min¬ 
imum around Tc and zero pb [47H5()j . The smaller the 
collision energy, the larger the average baryon chemical 
potential in the system. This indicates that the physical 
value of 77 /s should increase with increasing ps- 

In Ref. [5T] it was argued that 77 /s is not an ap¬ 
propriate measure of the fluidity of the system. How¬ 
ever, the measure of fluidity proposed in that paper, 
Ljj/Ln = /{wCs), where n is the total particle 

number density, w enthalpy, and the speed of sound, is 
difficult to implement in the present fluid-dynamical cal¬ 
culation since n is not well defined in our two-phase EoS. 
Instead, we use as an alternative measure of fluidity the 
combination pT/w = pT/{e + P) = p/{s + J2a fJ-ana/T), 
where Ua are the charge densities (baryon, strange, elec¬ 
tric) and Pa the corresponding chemical potentials, and 
which approaches p/s in the limit of small charge den¬ 
sities. We have performed an additional round of sim¬ 
ulations, keeping pT/w = 0.08 and p/s = 0.08 at all 
collision energies to see whether different measure of flu¬ 
idity makes any difference. The resulting elliptic and 
triangular flow coefficients are shown in Fig.|T3| One can 
see that at all considered collision energies there is no 
visible difference in the elliptic flow coefficient between 
the 77 /s = 0.08 and pT/w = 0.08 cases. We have also 
checked that the two scenarios result in virtually same pp 
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FIG. 13: (Color online) Pt integrated elliptic and triangular 
flow coefficients V 2 and V 3 as a function of collision energy. 
Solid red line represents the results from Fig. |10| obtained us¬ 
ing collision energy dependent 77 /s. Dashed blue and dotted 
green lines correspond to collision independent pTjw — 0.08 
and 77 /s = 0.08, respectively. In all three cases the other 
model parameters were taken to depend on the collision en¬ 
ergy as shown in Ta ble [II| T he experimental data is from the 
STAR collaboration |10l |40 |. 


spectra and dN/dy distributions. This indicates that the 
contribution from baryon/electric charge density to the 
entropy density does not induce strong enough baryon 
density dependence of the 77 /s ratio to affect the hydro- 
dynamic evolution. 


RHIC energy, = 7.7 — 200 GeV. After tuning the 
parameters, it was possible to reproduce the observed 
pseudorapidity and transverse momentum distributions 
of produced hadrons and their elliptic flow coefficients. 
The reproduction of the data requires a finite shear vis¬ 
cosity over entropy density ratio rj/s which depends on 
collision energy. This ratio was found to decrease from 
77 /s = 0.2 to 0.08 as collision energy increases from 
y'sNN = 7-7 to 39 GeV, and to stay at y/s = 0.08 for 
39 < ^/s < 200 GeV. Since the average baryochemical 
potential at midrapidity decreases with increasing col¬ 
lision energy, the required collision energy dependence 
of the effective 77 /s indicates that the physical ? 7 /s-ratio 
may depend on baryochemical potential, and that 77 /s 
increases with increasing p,B- It was also found that a 
constant and collision energy independent pT/w = 0.08 
and 77 /s = 0.08 in hydrodynamic phase yield quantita¬ 
tively similar results. This indicates that the /xsttb term 
in entropy density does not induce the baryon density 
dependence of 77 /s required to reproduce the data when 
rjT/w \s kept independent of collision energy. 

In addition we have explored the parameter depen¬ 
dence of the model results and generally found a < 10% 
variation of the results, when the individual parame¬ 
ters were varied by 10%. Of course, the proper evalu¬ 
ation of the effect of finite baryochemical potential on 
77 /s would require reproducing all the data using the 
same temperature and baryochemical potential depen¬ 
dent parametrization of rj/s at all energies and centrali¬ 
ties. This will be addressed in future studies. 


VI. SUMMARY AND OUTLOOK 

A hybrid model featuring a 3-1-1-dimensional viscous 
hydrodynamic phase with an explicit treatment of finite 
baryon and charge densities is introduced. The model 
employs a chiral model equation of state for the hydrody¬ 
namic stage. The initial and late non-equilibrium stages 
are modeled using the UrQMD hadron cascade on an 
event-by-event basis. 

This hybrid model was applied to describe the dynam¬ 
ics of relativistic heavy ion collisions at energies ranging 
from the lowest RHIC beam energy scan energy to full 
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